Strong-coupling expansions for the pure and disordered bose 

Hubbard model 



J. K. Freericks^ and H. Moment . 
( a ) Department of Physics, Georgetown University, Washington, DC 20057 
( 6 ) Theoretische Physik, ETH Honggerberg, CH-8093 Zurich, Switzerland 

(February 1, 2008) 

Abstract 

A strong-coupling expansion for the phase boundary of the (incompressible) 
Mott insulator is presented for the bose Hubbard model. Both the pure case 
and the disordered case are examined. Extrapolations of the series expansions 
provide results that are as accurate as the Monte Carlo simulations and agree 
with the exact solutions. The shape difference between Kosterlitz-Thouless 
critical behavior in one-dimension and power-law singularities in higher di- 
mensions arises naturally in this strong-coupling expansion. Bounded disor- 
der distributions produce a "first-order" kink to the Mott phase boundary in 
the thermodynamic limit because of the presence of Lifshitz's rare regions. 
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I. INTRODUCTION 

Strongly interacting bosonic systems have attracted a lot of recent interestllHi Physi- 
cal realizations include short correlation-length superconductors, granular superconductors, 
Josephson arrays, the dynamics of flux lattices in type II superconductors, and critical be- 
havior of 4 He in porous media. The bosonic systems are either tightly bound composites of 
fermions that act like effective bosonic particles with soft cores, or correspond to bosonic 
excitations that have repulsive interactions. For this reason, these systems are modeled by 
soft-core bosons which are described most simply by the bose Hubbard model. Various as- 
pects of this model were investigated analytically by mean-field theoryS'i, by renormalization 
group techniquesB! and by projection methodsi. The bose Hubbard model has also been 
studied with quantum Monte Carlo methods (QMC) by Batrouni et al.i in one dimension 
(1+1) and by Krauth and TrivediS, van Otterlo and Wagenblasti, and Batrouni et al.i in 
two dimensions (2+1). In this contribution, the Mott phase diagram is obtained from a 
strong-coupling expansion that has the correct dependence on spatial dimensionality, is as 
accurate as the QMC calculations, and agrees with the known exact solutions. Preliminary 
results for the pure case have already appeared^. 

The bose Hubbard model is the minimal model which contains the key physics of the 
strongly interacting bose systems — the competition between kinetic and potential energy 
effects. Its Hamiltonian is 

ij i % i 

where b{ is the boson annihilation operator at site i, tij is the hopping matrix element 
between the site i and j, is the local site energy, U is the strength of the on-site repulsion, 
and [x is the chemical potential. The hopping matrix is assumed to be a real symmetric 
matrix (t^ = tji) and the lattice is also assumed to be bipartite; i. e., the lattice may be 
separated into two sublattices (the A sublattice and the B sublattice) such that £y vanishes 
whenever i and j both belong to the same sublattice (in particular, this implies tu = 0). 
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The local site energy is a quenched random variable chosen from a distribution of site 
energies that is symmetric about zero and satisfies J2i e i — 0- The pure case corresponds to 
all site energies vanishing (e, = 0). 

The form of the zero temperature (T = 0) phase diagram can be understood by starting 
from the strong-coupling or "atomic" limitJMl. In this limit, the kinetic energy vanishes 
(tij = 0) and every site is occupied by a fixed number of bosons, n . In the pure case, the 
ground-state boson occupancy (no) is the same for each lattice site, and is chosen to minimize 
the on-site energy. If the chemical potential, /i = (no + 5)U, is parameterized in terms of the 
deviation, 8, from integer filling n , then the on-site energy is E(n ) = —5Un — ^Un (n + l), 
and the energy to add a boson onto a particular site satisfies E(no + l) — E(uq) = —5U. Thus 
for a nonzero 5, a finite amount of energy is required to move a particle through the lattice. 
The bosons are incompressible and localized, which produces a Mott insulator. For 5 = 0, the 
ground-state energies of the two different boson densities are degenerate [E{tiq) = E(n + 1)] 
and no energy is needed to add or extract a particle; i. e., the compressibility is finite and 
the system is a conductor. As the strength of the hopping matrix elements increases, the 
range of the chemical potential fi about which the system is incompressible decreases. The 
Mott-insulator phase will completely disappear at a critical value of the hopping matrix 
elements. Beyond this critical value of the system is a superfluid. 

In the disordered case, a Mott-insulating phase may or may not exist depending upon the 
strength of the disorder. The energy to add a boson onto site i becomes E(rio + l) — E(rio) = 
€i — SU, so that the system is compressible if a site % can be found which satisfies = SU. 
If the disorder is assumed to be symmetrically bounded about zero (|e^| < AU) then a 
Mott insulator exists whenever A < |. The ground-state boson occupancy is uniformly 
equal to n within the Mott insulating phase which extends from —A > 5 > A — 1 (when 
tij = 0). Once again, the bosons are incompressible within the Mott phase and the system is 
insulating. As the hopping matrix elements increase in magnitude, the range of the chemical 
potential within which the system is incompressible decreases until the Mott phase vanishes 
at a critical value of the hopping matrix elements. The compressible phase will typically 
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also be an insulator and is called a bose glassB, but it has been conjectured that in some 
cases the transition proceeds directly from the Mott insulator to the superfluid0i. 

The phase boundary between the incompressible phase (Mott insulator) and the com- 
pressible phase (superfluid or bose glass) is determined here in a strong-coupling expansion 
by calculating both the energy of the Mott insulator and of a defect state (which con- 
tains an extra hole or particle) in a perturbative expansion of the single-particle terms 
— J2ijtijb\bj + J2i e ini- At the point where the energy of the Mott state is degenerate with 
the defect state, the system becomes compressible. In the pure case, the compressible phase 
is also superfluid, but in the disordered case, the compressible phase is a bose glass (except 
possibly at the tip of the Mott lobe)0i. 

There are two distinct cases for the defect state: 8 < corresponds to adding a particle to 
the Mott-insulator phase (with hq bosons per site); and S > corresponds to adding a hole 
to the Mott-insulator phase (with no + 1) bosons per site. Of course, the phase boundary 
depends upon the number of bosons per site, n , of the Mott insulator phase. 

To zeroth order in tij/U the Mott-insulating state is given by 

N -i 

t=x vW v ' 

where no is the number of bosons on each site, N is the number of sites in the lattice and 
|0) is the vacuum state. The defect state is characterized by one additional particle (hole) 
which moves coherently throughout the lattice. To zeroth order in the single-particle terms 
the wave function for the "defect" state is determined by degenerate perturbation theory: 

l*DcfK)> h ° ie = 4=E i*Mott(n )> (0) (3) 

V n o i 

where the /j is the eigenvector of the corresponding single-particle matrix S^J (no) = —tij+ 
3ij£i/( n o + 1) [Sif (no) = —tij — Sijei/no] with the lowest eigenvalue (the hopping matrix 
is assumed to have a nondegenerate lowest eigenvalue). It is well known that the minimal 
eigenvalue of the single-particle matrix Sij is larger than the sum of the minimal eigenvalue 
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of the hopping matrix plus the minimal eigenvalue of the disorder matrix. However, it 
has been demonstrated that as the system size becomes larger and larger, the minimal 
eigenvalue approaches the sum of the minimal eigenvalues of the hopping matrix and of the 
disorder matrix as closely as desiredS (because of the existence of arbitrarily large "rare 
regions" where the system looks pure with e« = —AU or with = AU). Therefore, in the 
thermodynamic limit, the perturbative energy of each defect state becomes 

4 P e f Vo) " ^MottK) = -8^U + A min (n + 1) - AU + ... (4) 

4 h c°f l0) K) - ^Mott(no) = 5 {holc) U + \ mhl n - AU + ... (5) 

to first order in S, where A m j n is the minimal eigenvalue of the hopping matrix — t^. In 
the case of nearest-neighbor hopping on a hypercubic lattice in d- dimensions, the number 
of nearest neighbors satisfies z = 2d and the minimal eigenvalue is A m i n = —zt. 

The boundary between the incompressible phase and the compressible phase is deter- 
mined when the energy difference between the Mott insulator and the defect state vanishes 
(the compressibility is assumed to approach zero continuously at the phase boundary). Thus 
two branches of the Mott lobe can be found depending upon whether the defect state is an 
additional hole or an additional particle. The two branches of the Mott-phase boundary 
meet when 

5 (part) K) + l = 5 {holc \n ). (6) 

The additional one on the left hand side arises because S is measured from the point fi/U = 
uq. Equation may be used to estimate the critical value of the hopping matrix element 
beyond which no Mott-insulator phase exists. Let x denote the combination dt/U and 
consider the first-order expansions in Eqs. (f|) and (||). The critical value of x satisfies 

which vanishes when the disorder strength becomes too large (A > 1/2). Note that the 
critical value of x is independent of the dimension of the lattice; the dimensionality first 



enters at second order in t. The slope of the phase boundaries about the point \i = tiqU are 
equal in magnitude [lim^o ^<5( part )(no, x) = —\im x ^ -^S^ hole \no + l,x)], but change their 
magnitude as a function of the density n , implying that the Mott-phase lobes always have an 
asymmetrical shape. Note further that the presence of disorder shifts the phase boundaries 
uniformly by A, but the slope is independent of the disorder distribution. 

The bose Hubbard model in the absence of disorder is examined by a strong-coupling 
expansion through third order in the single-particle matrix S in Section II. The exact solution 
for an infinite-dimensional latticJ is examined and various different extrapolation techniques 
are employed that do and do not utilize additional information of the scaling analysis of the 
critical point. Section III describes the similar results for the disordered bose Hubbard model 
and a discussion follows in Section IV. 



II. THE PURE CASE 



The bose Hubbard model in Eq. (|1|) is studied in the absence of disorder (ej = 0). The 
many-body version of Rayleigh-Schrodinger perturbation theory is employed throughout. 
To third order in tij/U, the energy of the Mott state with no bosons per site becomes 



^Mott(no) = N 



1 1 t 2 - 

-5Un - -Un (n + 1) - —J2 jj n o( n o + 1) 

ij 



(8) 



which is proportional to the number of lattice sites N. Note that the odd-order terms in 
tij/U vanish in the above expansion (odd-order terms may enter for nonbipartite lattices). 
The energy difference between the Mott insulator and the defect state with an additional 
particle (5 < 0) satisfies 

4 P ef rt) (no) - E Mott (n ) = -S&^U + A min K + 1) + 217 Eii t 2 Jfn (5n + 4) - ^A^oK + 



+ 772^0(^0 + 1) 



(2n + l)X 3 mhl - (f n + |) A min £y t%f] - (An + 2) f,^]) 



to third order in tij/U; while the energy difference between the Mott insulating phase and 
the defect phase with an additional hole (5 > 0) satisfies 
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£&°f le) (no) - ^Mott(no) = 5 (holc) f7 + A min n + ^ £ y + l)(5n + 1) - ^A^oK + 1) 



+ 772^0(710 + 1) 



(2n + 1)A 
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min 



f n + ^)A min Eij *?•/? - (4n + 2) ^ 



The eigenvector is the minimal eigenvector of the hopping matrix — tij with eigenvalue 
A m i n and is identical in the particle and hole sectors. These results have been verified by 
small-cluster calculations on two and four-site clusters. Note that the energy difference in 
Eqs. ([|) and fllPf ) is independent of the lattice size iV indicating that QMC simulations 
should not have a very strong dependence on the lattice size. 

In the case of nearest-neighbor hopping on a d- dimensional hypercubic lattice, the mini- 
mum eigenvalue satisfies A min = —zt, the sum J2ij tfjfj becomes zt 2 , and the sum Ey Mijfj 



equals zt 3 . Equations (^[) and (|I0D can then be solved for the shift in the chemical potential 
5 at which the system becomes compressible as a function of the parameter x = dt/U. The 
results for the upper boundary are given by 



ii , x) = -2x(n + 1) + ^x 2 n (5n + 4) - Ax 2 n (n + 1) 



+ 2x 3 n (n + 1) 



:-s + 1 - JrH + 



-4-1-2 — 



(11) 



to third order in x, and the lower boundary is given by 



6^ hole \n ,x) 



2xn — ^x 2 (n + l)(5n + 1) + Ax 2 n (n + 1) 



- 2x 3 n (n + T 



+ £-£)«o + (-4 + &- 



2d 
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2d d 2 > 



(12) 



to third order in x. 

As a further check on the accuracy of the Mott phase boundaries in Eqs. (|Tl|) and fljj), 
we compare the perturbative expansion to the exact solution on an infinite-dimensional 
hypercubic latticeB (which corresponds to the mean-field solution). Note that the solution 
in Ref. [l| was for the infinite-range-hopping model; this solution is identical to that on an 
infinite-dimensional lattice in the pure case. The Mott phase boundary may be expressed 
as 



no 

U 



- x ± J x 2 - x(2n + 1) + - 



(13) 



where the plus sign denotes the upper branch to the Mott lobe (<5( part )), and the minus sign 
corresponds to the lower branch (£( hole ) — 1). The critical point can also be determined as 
the value of x where the square root vanishes. One finds 



which depends on n Q as l/n in the limit of large n Q . The strong-coupling expansions flllf) 
and (|l^) agree with the exact solution fll3|) when the latter is expanded out to third order 
in x, providing an independent check of the algebra. Note further, that the exact solution 
uniquely determines the expansion coefficients of the powers of x that do not involve inverse 
powers of d and the perturbation expansion is only required to determine the 1/d corrections. 

The strong-coupling expansion for the x, \x phase diagram in one dimension is compared 
to the QMC results of Batrouni et al.i in Figure 1. The solid lines indicate the phase 
boundary between the Mott-insulator phase and the superfluid phase at zero temperature 
as calculated from Eq. ([11]) and Eq. (|T2|). The solid circles are the results of the QMC 
calculation!! at a small but finite temperature (T ~ U/2). The dotted line is an extrapolation 
from the series calculation that will be described below. Note that the overall agreement 
of the two calculations is excellent. For example, the critical value of the hopping matrix 
element for the first Mott lobe (no) is x cr j t = 0.215, while the QMC calculations foundi 
^crit = 0.215 ± 0.01. A closer examination of Fig. 1 shows a systematic deviation of the 
lower branch for larger values of x. We believe that this is most likely a finite-temperature 
effect, since the Mott-insulator phase becomes more stable at higher temperatures!, and the 
systematic errors of the QMC calculation due to finite lattice size and finite Trotter error 
are easily controlled^. 

It is known from the scaling theory of Fisher et all that the phase transition at the tip of 
the Mott lobe is in the universality class of the (d + 1) dimensional XY model. Although a 
finite-order perturbation theory cannot describe the physics of the tricritical point correctly, 
we find that the density fluctuations dominate the physics of the phase transition even close 
to the tricritical point. Note how the Mott lobes have a cusp-like structure in one dimension, 




(14) 
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mimicking the Kosterlitz-Thouless behavior of the critical point. 

Figure 2 (a) presents the strong-coupling expansion for the x, \x phase diagram in two 
dimensions. For comparison, the tricritical point of the first Mott-insulator lobe as obtained 
by the QMC simulations of Krauth and Trivedi@ is marked by a solid circle with error 
bars (the chemical potential for the tip of the Mott lobe was not reported in Ref. |7], so 
we fixed it to be /i C r«t)- The solid line is the strong-coupling expansion truncated to third 
order, while the dotted line is an extrapolation described below. Their simulation gives a 
critical value of x cr i t = 0.122 ± 0.006, whereas our calculation yields x cr i t ~ 0.136 which is 
in reasonable agreement. Note that the qualitative shape of the Mott lobes has changed 
from one dimension to two dimensions, mimicking the power-law critical behavior of the 
XY model in three or larger dimensions. 

Figure 2 (b) shows the corresponding figure for the Uq — > oo limit corresponding to 
the quantum rotor model. The QMC results are from van Otterlo and Wagenblasti. The 
horizontal axis has been rescaled to = lim n0 _ ) . oo UqX. We believe that the relatively large 
difference between the QMC and the strong-coupling perturbation theory arises from the 
use of the Villain approximation in the QMC simulations. 

Finally the strong-coupling expansion is compared to the exact calculation in infinite 
dimensions^. In infinite dimensions, the hopping matrix element must scale inversely with 
the dimensionlll, t = t*/d, t* = finite, producing the mean-field-theory result of Eq. fll3"l). In 
Figure 3 the strong-coupling expansion (solid line) is compared to the exact solution (dashed 
line) and to an extrapolated solution (dotted line) which will be described below. Even in 
infinite dimensions, the agreement of the strong-coupling expansion with the exact results 
is quite good. 

As a general rule, the truncated strong-coupling expansions appear to be more accurate 
in lower dimensions, which implies that the density fluctuations of the bose Hubbard model 
are also more important in lower dimensions. 

At this point we turn our attention to techniques which enable us to extrapolate the 
strong-coupling expansions to infinite order in hopes of determining a more accurate phase 



diagram. The simplest method is called critical-point extrapolation. The critical point 
(yU cr it, Xcnt) is calculated at each order (m) of the strong-coupling expansion and is extrapo- 
lated to infinite order (m — > oo). The ansatz that the extrapolation is linear in 1/m can be 
checked by determining the correlation coefficient r of the critical points (a value of \r\ that 
is near 1 indicates a linear correlation). The correlations are found to be most linear for large 
dimensions (|r| = 0.99999 in infinite dimensions for the first Mott lobe) but remain fairly 
linear even in one dimension (|r| > 0.995 for the x cr i t extrapolation and |r| > 0.95 for the 
/x cr it extrapolation). Since the second- and third-order expansions are expected to be more 
accurate than the first-order calculation, we adopt the following strategy for performing the 
extrapolations: the results of the second- and third-order expansions are extrapolated to 
m — > oo to determine the estimate for the critical point, and the results of the first, second 
and third orders are then extrapolated to m — > oo in order to estimate the error in the 
critical point, and to test the linear-extrapolation hypothesis. The error estimate is chosen 
to be 1.5 times as large as the difference between the two different extrapolations. Figure 4 
plots the critical hopping matrix elements x cr ;t versus 1/m for the infinite-dimensional case 
and n = 1,2,3. The solid dots are the results of the strong-coupling expansion truncated 
to mth order and the solid line is the linear extrapolant. The open circles are the exact 
solutions from Eq. flli]). Note that although the linear correlation coefficient is very close to 
1, the error in the critical point is about 2%. The results for the critical-point extrapolation 
are recorded in Table I. 

The critical-point extrapolation does not yield any information on the shape of the Mott 
lobes, but only determines the critical point. An alternate extrapolation technique, called 
the chemical-potential extrapolation method will determine an extrapolated Mott-phase 
lobe and critical point. The idea is to fix the magnitude of the hopping matrix elements 
and determine the value of the chemical potential from Eqs. (|TTD and ( |T2"D for the upper 
and lower branch of the Mott lobe. The chemical potential is determined from a first, 
second, and third-order calculation and then extrapolated to infinite order assuming the 
ansatz of a linear dependence upon 1/m. This procedure determines an extrapolated Mott 
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lobe that should be more accurate than the truncated strong-coupling series. The result 
for the infinite-dimensional case is presented as a dotted line in Fig. 3. Note that the 
critical point is not determined as accurately by this technique as it was in the critical-point 
extrapolation method. The chemical-potential extrapolation method fails in one dimension 
since the extrapolated branches of the extrapolated Mott lobe do not close. 

A third approach is to use the results of the scaling theory^. The critical point is that 
of a (d + l)-dimensional XY model, and therefore, has a Kosterlitz-Thouless shape in one 
dimension and a power-law shape in higher dimensions. Examination of the exact result for 
infinite dimensions flUf ) leads one to propose the following ansatz for the Mott lobe in d > 2 

jj-n = A(x) ± B{x)(x CTit - xT (15) 

where A(x) = a + bx + cx 2 + ... and B(x) = a + (3x + r yx 2 + ... are regular functions of x (that 
should be accurately approximated by their power-series expansions) and zv is the critical 
exponent for the (d + l)-dimensional XY model. In the unconstrained-scaling-analysis 
extrapolation method the exponent zv is determined by the strong-coupling expansion in 
addition to the parameters a, b, c and a, (3, 7. This provides a perturbative estimate of 
the exponent zv which can be checked against its well-known values. In the constrained- 
scaling-analysis extrapolation method zv is fixed at its predicted valuesB of zv ~ 2/3 in two 
dimensions and zv = 0.5 in higher dimensions. In direct analogy to Eq. (|i~5l) , we propose 
the Kosterlitz-Thouless form 



jj — n = A(x) ± B(x) exp 



W 



(16) 



V ^crit 

for the constrained-extrapolation-method in one dimension. 

When the unconstrained-scaling-analysis extrapolation method is carried out, one finds 
that there is no solution for the critical exponent in one dimension (which is consistent with 
Kosterlitz-Thouless behavior), that in d = 2 the exponent satisfies zv w 0.58, in d = 3 the 
exponent is zv w 0.54, and in infinite dimensions zv ~ 0.5. There is a slight no dependence 
to the exponents that are calculated in this method, but that arises from the truncation of 
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the series to such a low order. In general, the unconstrained extrapolation method produces 
an accuracy of about 15% in the exponent zv, and the method appears to work best in 
higher dimensions. 

The results for the constrained-extrapolation method are plotted with a dotted line in 
Fig. 1 for the one- dimensional case. The values of the critical points are (^ cr u = 0.186, x cr u = 
0.265), (1.319,0.155), and (2.371,0.111) for n = 1,2,3, respectively. These critical points 
occur at larger values of x than predicted by the QMC simulations!, but it is difficult 
to gauge whether the extrapolated series expansions are more or less accurate than the 
QMC simulations because of the finite-temperature effects in the latter. The constrained- 
extrapolation results in two dimensions are plotted with a dotted line in Fig. 2 (a) and 
(b). The values of the critical points are (0.375,0.117), (1.426,0.069), (2.448,0.049) for 
no = 1, 2, 3, respectively. The agreement with the QMC simulationsB is excellent. Similarly, 
the extrapolated critical point for the 2 — d rotor model is (0.5,0.171) which also agrees 
well with the QMC. In this latter case [Fig. 2 (b)] the errors between the extrapolated 
series expansion and the QMC can be traced to the use of the Villain approximation in 
the latter. When the constrained-scaling-analysis extrapolation method is applied to the 
infinite-dimensional case the result is indistinguishable from the exact solution when the 
two are plotted on the same graph. 

The extrapolation techniques work best in higher dimensions virtually producing the 
exact result in infinite-dimensions. This gives us confidence that the extrapolated results 
of the series expansions can produce numerical answers that are at least as accurate as the 
QMC simulations. 

III. THE DISORDERED CASE 

The most common type of disorder distribution that has been considered in relationship 
to the "dirty" boson problem is the Anderson model (continuous disorder distribution). In 
the Anderson model the distribution p(e) for the on-site energies {ej} is continuous and flat, 
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satisfying 



p(e) =9(A-e)9(A + e) 



1 



(17) 



2A 



The symbol A denotes the maximum absolute value that the site energy can assume for a 
given (bounded) distribution (|ej| < AU). This disorder distribution is symmetric [p(— e) = 
p(e)] and in particular it satisfies J2i e i — 0. The results presented in this contribution are 
insensitive to the actual shape of the disorder distribution; all we require is a symmetric 
distribution with |ej| < AU. 

We begin by reexamining the exact solution of the infinite-range-hopping mo delB. If all 
energies are measured in units of the boson-boson repulsion U, then the analysis of Ref. [I] 
derives an equation that relates the hopping matrix elements to the chemical potential at 
the Mott phase boundary 



with y = —uq + h/U the chemical potential and p(e) the disorder distribution. This solution 
assumes that the phase transition from the Mott phase to the bose glass is a second- order 
phase transition. 

When Eq. (|i~8|) is solved for the Anderson model distribution fllTD , one finds that 
the chemical potential for the lower branch of the Mott lobe behaves like y = —A — 
exp[— l/2x(no + 1)] for small x. This result is nonperturbative in the hopping matrix el- 
ements and cannot be represented by a simple perturbative theory about x = 0. The reason 
why this happens is that the infinite-range-hopping model has no localized states for any 
disorder distribution (however, this statement does depend on the disorder distributional) . 
Localized states can occur in the infinite-dimensional limit at the tails of the distribution. 
Therefore we expect that the transition will have a different qualitative character on a hyper- 
cubic lattice with nearest-neighbor interactions. In fact, the perturbative arguments given 
in the introduction, show that the phase boundaries have the same slope as they did in 
the pure case. Furthermore, we expect that the transition to be first order at the tip of 




(18) 
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the Mott lobe because the states that the bosons initially occupy in the compressible phase 
are localized within the rare regions of the lattice (where the site energies are all equal to 
—AU) implying that there is no diverging length scale at the transition. The perturbative 
expansion for the energy of the Mott phase is unchanged from Eq. (|8|) in the presence of 
disorder (if the disorder distribution satisfies — 0), while the defect phases have a 

trivial dependence upon disorder (in the thermodynamic limit) — the energy for a particle 
or hole defect state is shifted by — AU, so the effect of the disorder is simply to shift the 
Mott-phase boundaries inward by AU. The critical point, where the Mott phase disappears, 
is no longer described by a second-order critical point [in which the slope of fi(x) becomes 
infinite at x cr n\ but rather is described by a first-order critical point [in which the slope of 
fj,(x) changes discontinuously at x cr i t }. 

In the thermodynamic limit one can always find a rare region of arbitrarily large extent 
which guarantees the existence of the first-order transition, but the density of these rare 
regions is an exponentially small function of their size. For this reason the compressibility 
at the Mott-phase boundary will also be exponentially small. Finite-size effects play a much 
more important role in the disordered case: it is impossible to determine the Mott-phase 
boundary accurately in the thermodynamic limit by scaling calculations performed on small 
lattices, because the lattice size must be large enough to contain rare regions within which 
the bosons can be delocalized. (Finite-size effects can be studied with the strong-coupling 
expansion which is given to third order in the single-particle matrix S in the appendix.) 

The most accurate way of calculating the Mott phase boundary is then to take the results 
of the constrained-scaling-analysis extrapolation for the pure case and shift the branches by 
the strength of the disorder. This is plotted in Figure 5 for the one-dimensional case and 
two different values of disorder (A = 0,0.25). The thermodynamic limit is represented by 
the solid line for the pure case, and dotted lines for the disordered case, while the dashed 
line is the result of an Anderson-model disorder distribution on a finite lattice with 256 sites. 
The QMC results of Batrouni, et. a/.i correspond to lattice sizes ranging from 16 sites to 
256 sites (the disorder parameter is A = for the solid dots and A = 0.25 for the open 
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dots). The Mott phase is stabilized on finite-sized systems because they do not possess the 
rare regions needed to correctly determine the Mott phase boundary This is clearly seen in 
the QMC results, which predict a much larger region for the Mott phase than the strong- 
coupling perturbation theory does in the thermodynamic limit. The perturbative results for 
a finite system are much closer to the QMC results as expected. (Note that the finite-size 
calculations have not been extrapolated, so they should underestimate the stability of the 
Mott phase in one dimension which is clearly seen in Figure 5.) Also the slope of the phase 
boundary approaches zero (as x — > 0) for the finite-size systems^. 

Figure 6 plots the Mott-phase diagram for the disordered bose Hubbard model in two 
dimensions and one value of the disorder (A = 0,0.182) in the thermodynamic limit. The 
solid dot (with error bars) is the QMC resultB (for the pure case with A = 0) and the 
open dot is the disordered case (A = 0.182). Note that in dimensions larger than one, the 
finite-size effects for the tip of the Mott lobe are not as strong as they are in one dimension. 

We compare in Figure 7 the differences between the infinite-dimensional lattice and the 
infinite-range-hopping model of Ref. [TJ The solid lines correspond to the exact solution with 
no disorder, the dotted lines are the infinite-dimensional lattice strong-coupling expansion 
with disorder (A = 0.2), and the dashed lines are the exact solution of the infinite-range- 
hopping model. In the case of disorder, the first-order nature of the transition is evidenced 
by the jump in the slope of the Mott phase boundary at the tip of the lobe. The second-order 
phase boundaries predict a more stable Mott phase, and their slopes all approach zero as x — > 
0. We expect in the region in between the (first-order) infinite-dimensional phase boundary 
and that of the (second-order) infinite-range-hopping model that the compressibility will be 
exponentially small, and will only become sizable as the second-order phase boundary is 
crossed. 

Because the Mott-phase to bose-glass phase transition is first order for the disordered 
case, and since the bose-glass to superfluid transition is always second order (because it 
involves a collective excitation that extends through the entire lattice), it is quite unlikely 
that there would ever be a region where the Mott phase has a transition directly to the 
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superfluid. The presence of the Lifshitz rare regions strongly supports the picture that the 
Mott phase is entirely enclosed within the bose-glass phase. This result is independent of 
any perturbative arguments, since the rare regions must dominate the Mott to bose-glass 
transition in the exact solution too. 

Finally, we calculate the dependence of the critical value of x at the tip of the Mott-phase 
lobe, as a function of the disorder strength A. Figure 8 plots this value of x cr i t (A)/x cr it(0) 
versus A for the one-, two-, and infinite-dimensional cases. The plot is limited to the lowest 
Mott-phase lobe with n — 1. Since the one-dimensional Mott phase lobes have a cusp- 
like behavior that is removed when disorder is added to the system, we expect x cr u to 
decrease very rapidly for small disorder. This effect is sharply reduced in higher dimensions. 
In the strong-disorder limit, the phase diagram is dominated by the first-order terms in 
the perturbative expansion, which have a trivial dependence on the dimensionality, but 
the slopes curves of the curves are unequal because a? C rit(0) depends strongly upon the 
dimensionality. 

IV. CONCLUSION 

We have developed a strong-coupling perturbation-theory approximation to the bose 
Hubbard model on a bipartite lattice. The perturbative results can be extrapolated in a 
number of different ways which either do or do not take into account the scaling theory of 
the critical point at the Mott tip. We find that a perturbative expansion thru third order 
rivals the accuracy of the QMC simulations, and is likely the best method for determining 
the Mott phase boundary of these interacting bose systems. 

We treated two different cases: the pure case and the disordered case. In the pure case 
the tip of the Mott lobe satisfies a scaling relation because it corresponds to a second-order 
phase transition in a d + 1-dimensional XY model. This is because the compressible phase 
is also superfluid which implies there is a diverging length scale at the phase transition. 
Calculations in the pure case are insensitive to finite-size effects. In the disordered case 
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we argued that the tip of the Mott phase lobe corresopnds to a first-order phase transition 
because the initial single-particle excitations are localized into the rare regions of the Lifshitz 
tails for any bounded disorder distribution. As a result there is a kink in the Mott phase 
boundary since the slope of /i(x) has a discontinuous jump at the tip of the lobe. In this case, 
there are strong finite-size effects because "typical" disorder distributions on finite lattices 
do not have Lifshitz tails. 

The results of these perturbative calculations have been compared to the available QMC 
simulations. We find a remarkable agreement between the two and are unable to determine 
which method is more accurate in a quantitative determination of the phase boundaries. 

The perturbation theory described here falls short in one aspect — it is unable to deter- 
mine the bose-glass-superfluid phase transition in the disordered case. It is possible that 
such a calculation could be performed, but since the particle density at which it occurs is 
a priori not known, and since the superfluid susceptibility diverges in the bose-glass phase, 
such a calculation may be problematic. 
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APPENDIX A: FINITE-SIZE EFFECTS FOR THE DISORDERED CASE 



The Rayleigh-Schrodinger perturbation theory for the case with disorder is straightfor- 
ward, but rather tedious. In the thermodynamic limit, the perturbative expansion simplifies 
because it is dominated by the rare regions of the lattice. For a finite-sized system, the 
perturbative results are more complicated. We summarize here the main results for a third- 
order strong-coupling expansion in the presence of disorder. 

The only assumptions that are made about the disorder distribution is that it is bounded 
and symmetric, so that |e«| < AU and Ei Q = 0. Restriction is also made to bipartite lattices. 

The upper phase boundary for the Mott to bose-glass transition is found by solving the 
equation 

«j(p° rt )(n ) = £H + 1) + !§^ (5n + 4) - E ijk /.^Mo(«o + 1) 



- n (n + 1 

Zt 2 ft 



E ijk i ti-^mno + 2) + A E . jk f^f k ( no + 1) 
W Zij fifWrio + 6) + £ Eij /i^/i(4no + 2) - A^(> + |) 

"o £y ||/|(I^o + 2) - f£ % ft Cino + 1) 



( e i~ e j+ e k) £ Ujtjk £ 

Ajk u J i jj2 J k 
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where fa is the minimal eigenvector of the single-particle matrix = —tij + €i5ij/ (n + 1) 

and A is its corresponding eigenvalue. The identity Ei? tijfj — Zt 2 was needed in deriving 
the above result. A similar calculation yields 

5( hole \n ) = -A no - \%{n Q + l)(5n + 1) + E^^S^oK + 1) 
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for the lower phase boundary of the Mott to bose-glass transition. Here we have that is 



the minimal eigenvector of the single-particle matrix 
corresponding eigenvalue. 



(hole) 



-tij — €i5ij/no and A is its 



In the thermodynamic limit we know that the minimal eigenvalue occurs in the rare 
regions where the disorder is constant and equal to its extreme value. The ground-state 
eigenvector is delocalized within the rare region (to minimize it's kinetic energy) and lo- 
calized to the rare region (to minimize its disorder energy). Such an eigenvector is now 
separately an eigenvector of the kinetic-energy matrix and of the disorder matrix, so we have 
— YlijUjfj = A/j, J2j e i$ijfj — —AUfi, and A = A — AU/(tiq + 1) with similar formulae 
for the g eigenvector. Plugging these thermodynamic limits into Eqs. ( |Al| ) and (|A2|) then 
yields the result that the Mott phase boundaries are only shifted uniformly by AU, and all 
higher-order dependence on the disorder vanishes. 
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TABLES 

TABLE I. Results for the critical-point extrapolation method described in the text. The 
critical point is recorded for the first three Mott lobes in one, two, three and infinite dimensions. 
Where possible the results from other calculation techniques are displayed in the last column. 



dimension 


n 


/^crit 


•^crit 


x crit (exact) 


1 


1 


0.255±0.11 


0.245±0.012 


0.215±0.010 a 




2 


1.359±0.06 


0.145±0.009 


0.130±0.020 a 




3 


2.400±0.04 


0.103±0.006 


0.104±0.020 a 


2 


1 


0.388±0.05 


0.114±0.013 


0.122±0.006 b 




2 


1.435L0.03 


0.067±0.006 






3 


2.454±0.02 


0.048±0.004 




3 


1 


0.400±0.03 


0.101±0.008 






2 


1.441±0.02 


0.060±0.004 






3 


2.458±0.01 


0.042±0.002 




oo 


1 


0.416±0.001 


0.0843±0.001 


0.0858 c 




2 


1.451±0.002 


0.0494±0.002 


0.0505 c 




3 


2.465±0.001 


0.0351±0.001 


0.0359 c 



a Ref. 2 (Monte Carlo simulation at finite temperature) 
b Ref. 7 (Monte Carlo simulation at finite temperature) 
c Ref. 1 [Exact solution from Eq. (14)] 
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FIGURES 

FIG. 1. The x, jj, phase diagram of the bose Hubbard model in one dimension [d = 1). 
The solid lines give the phase boundaries of the Mott insulator to the superfluid state as deter- 
mined from a third-order strong-coupling calculation. The dotted line is the constrained fit to a 
Kosterlitz-Thouless form. The circles are the result of the QMC calculation of Batrouni et. al. 2 . 

FIG. 2. (a) The x, fi phase diagram of the bose Hubbard model in two dimensions (d = 2). 
The solid lines give the phase boundaries of the Mott insulator to the superfluid state as determined 
from a third-order strong-coupling calculation. The dotted line is the constrained fit to a power 
law form with exponent zv = 2/3. The point (with error bars) indicates the tricritical point as 
determined by the QMC calculation of Krauth and Trivedi 7 (no value for /j, cr i t was given in Ref. 7) . 
(b) The yco, fx phase diagram for the quantum rotor model in two dimensions. The solid lines are 
the perturbative results to third order and the dotted lines are the constrained extrapolation fit. 
The dots are the QMC results of van Otterlo and Wagenblast 8 . The disagreement between the 
QMC and the extrapolated results most likely arise from the use of the Villain approximation in 
the former. 

FIG. 3. The x, u phase diagram of the bose Hubbard model in infinite dimensions (d — > oo). 
The solid lines give the phase boundaries of the Mott insulator to the superfluid state as determined 
from a third-order strong-coupling calculation. The dashed lines are the result of the mean-field 
calculation of Fisher et. al. 1 . The dotted lines are the chemical-potential extrapolation described 
in the text. 

FIG. 4. The critical-point extrapolation method in infinite dimensions. The solid circles 
are the results for x cr it calculated at first, second, and third order. The solid line is the linear 
extrapolation to infinite order. The open circles are the exact result. Three cases are shown 
no = 1 (top), no = 2 (middle), and no = 3 (bottom). Note that although the three points have a 
correlation coefficient larger than 0.9999 the accuracy of the extrapolated critical point is only on 
the order of 2%. 
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FIG. 5. Phase diagram in one dimension with disorder. The perturbative approximations are 
plotted with solid (dotted) lines for the pure (disordered) cases in the thermodynamic limit. Three 
cases of disorder are included: (A = 0.125, 0.25, 0.375). The dashed line is the perturbative results 
for a finite system with 256 lattice sites. The solid dots are the quantum Monte Carlo results from 
Ref. 2 in the pure case, and the open dots are for the disordered case with A = 0.25 (the disordered 
calculations were performed on lattice sizes ranging from 16 to 256 sites). 

FIG. 6. Phase diagram in two dimensions with disorder. The perturbative approximations 
are plotted with solid (dotted) lines for the pure (disordered) cases in the thermodynamic limit. 
The disorder was set equal to A = 0.182. The solid (open) dots are the quantum Monte Carlo 
results of Ref. 7 for the pure (disordered) cases. Surprisingly, in two and higher dimensions, the 
finite-size effects in the disordered regime appear to be weaker. 

FIG. 7. Phase diagram in infinite dimensions with disorder. The perturbative approximations 
are plotted with solid (dotted) lines for the pure (disordered) cases in the thermodynamic limit. The 
disorder was set equal to A = 0.2. The dashed line is the exact solution of the infinite-range-hopping 
model from Ref. 1. Note how the Mott phase is more stable with the infinite-range calculation, 
and how it has vanishing slope as x — > 0. Interestingly, the location of the tip of the Mott lobe is 
close for both the infinite-dimensional and the infinite-range calculations. 

FIG. 8. Plot of x cr it(A)/x cr it(0) as a function of dimension. The solid line is for one dimension, 
the dotted line for two dimensions, and the dashed line for infinite dimensions. Note that the intitial 
decrease of x cr u is very rapid in one dimension, because of the cusp-like shape of the Mott lobe in 
the pure case, but is much slower in higher dimensions (because the tip has a power law behavior 
in the pure case). 
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